'''
# encoding=utf-8
@version: v1.0
@author: Jack Jiang
@contact: 2276009941@qq.com
@site: https://www.ms.com/
@software: PyCharm Community Edition
@file: 奇异值分解.py
@time: 2019-03-14 20:04
'''
import warnings

warnings.filterwarnings('ignore')


import numpy as np
from numpy.linalg import svd,eig

A = np.array([[126,  52,  -3, -69],
       [ 52, 292, -73, -80],
       [ -3, -73, 141, -31],
       [-69, -80, -31,  78],
       [-69, -80, -31,  178]])

print('***************A矩阵，行数m=5,列数n=4*************************')

# 奇异值分解: A.shape=(m,n) U.shape=(m,m) V.shape(n,n) S.shape=(m,n)
U, S, VT = svd(A)
V = VT.T
print('U=\n',U)       #左奇异矩阵U, U.shape=(m,m)
print('V=\n',V)    #右奇异矩阵V, V.shape=(n,n)

print('S=\n',S)    #奇异值
S = np.diag(S)     #转换成对角矩阵
print(S)
h = np.array([[0,0,0,0]])  # 补充一行0到对角矩阵
S = np.vstack((S,h))    #转换成奇异值矩阵S, S.shape=(m,n)
print(S)
print('A=\n',np.dot(np.dot(U, S), VT))    # A(m,n)  U(m,m)  Lambda(m,n)  V(n,n)

#U为单位正交矩阵
print('U模长=',((U**2).sum(axis=0))**0.5) #验证模长=>模为1
print('U向量内积=',np.dot(U[:,0],U[:,4])) #验证正交=>向量相互正交为0

#V为单位正交矩阵
print('V模长=',((V**2).sum(axis=0))**0.5) #验证模长=>模为1
print('V向量内积=',np.dot(V[:,0],V[:,2])) #验证正交=>向量相互正交为0

print('***************A矩阵，行数m=4,列数n=5*************************')

A=A.T  #行数m=4，列数n=5

# 奇异值分解: A.shape=(m,n) U.shape=(m,m) V.shape(n,n) S.shape=(m,n)
U, S, VT = svd(A)
V = VT.T
print('U=\n',U)    #左奇异矩阵U, U.shape=(m,m)
print('V=\n',V)    #右奇异矩阵V, V.shape=(n,n)

print('S=\n',S)    #奇异值
S = np.diag(S)     #转换成对角矩阵
h = np.array([[0],[0],[0],[0]])  # 补充一行0到对角矩阵
S = np.hstack((S,h))    #转换成奇异值矩阵S, S.shape=(m,n)
print('A=\n',np.dot(np.dot(U, S), VT))    # A(m,n)  U(m,m)  Lambda(m,n)  V(n,n)

#U为单位正交矩阵
print('U模长=',((U**2).sum(axis=0))**0.5) #验证模长=>模为1
print('U向量内积=',np.dot(U[:,0],U[:,3])) #验证正交=>向量相互正交为0

#V为单位正交矩阵
print('V模长=',((V**2).sum(axis=0))**0.5) #验证模长=>模为1
print('V向量内积=',np.dot(V[:,0],V[:,2])) #验证正交=>向量相互正交为0


print('证明右奇异矩阵V就是(A.T)乘以A后的矩阵特征值向量矩阵转置')
A=np.dot(A.T,A)
vals,vecs=np.linalg.eig(A)
print('vecs=\n',vecs)



# print('***************奇异值求解过程*************************')
# #验证奇异值分解过程：
# # (1) 先以A.T.dot(A)矩阵，求解特征值vals1和特征向量vecs1
# # (2) S=降序后的vals1的平方根)
# # (3) V=vecs1.T  (vecs1=vals1降序后对应的vecs1)
# # (4) 再以A.dot(A.T)矩阵，求解特征值vals2和特征向量vecs2
# # (5) U=vecs2    (vecs2=vals2降序后对应的vecs2)
#
# # (1) 先以A.T.dot(A)矩阵，求解特征值vals1和特征向量vecs1
# vals1,vecs1=eig(np.dot(A.T,A))
# # (2) S=降序后的vals1的平方根)
# print('S=sqrt(vals1)\n',np.sqrt(vals1))
# # (3) V=vecs1.T  (vecs1=vals1降序后对应的vecs1)
# print('V=vecs1.T\n',vecs1.T)
# # (4) 再以A.dot(A.T)矩阵，求解特征值vals2和特征向量vecs2
# vals2,vecs2=eig(np.dot(A,A.T));
# print(np.sqrt(vals2))  #vals1和vals2非零元素是相同的
# # (5) U=vecs2    (vecs2=vals2降序后对应的vecs2)
# print('U=vecs2\n',vecs2)
# # 将S转化为对角矩阵
# S=np.vstack((np.diag(np.sqrt(vals1)),np.array([[0,0,0,0]])))
# # 验证A=USV
# print('A=\n',vecs2.dot(S).dot(vecs1.T))